library(tidyverse)
library(readstata13)
library(foreign)
library(readxl)
library(writexl)
library(lubridate)
library(PanelMatch)

setwd("/Volumes/Backup Plus/Storage-OrganizedFiles/Article1-ReplicationMaterials_FullReplication/4.Diagnostics")
getwd()

load("OMC-ResultsMainModel_Updated.RData")

##############################################################################
### COBA 1
##############################################################################

##############################################################################
# Protest
##############################################################################

pdf(file = "Figure9.1-P.pdf",   
    width = 13, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 6))

coba_pm.none5 <- get_covariate_balance(pm.maha5b.p$att,
                                       data = data_final,
                                       use.equal.weights = TRUE,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       plot = TRUE,
                                       legend = F,
                                       ylim = c(-1, 1))

coba_pm.maha5b.p <- get_covariate_balance(pm.maha5b.p$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

coba_pm.ps.m5.p <- get_covariate_balance(pm.ps.m5.p$att,
                                         data = data_final,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.ps.w5.p <- get_covariate_balance(pm.ps.w5.p$att,
                                         data = data_final,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.CBPSm5.p <- get_covariate_balance(pm.CBPSm5.p$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

coba_pm.CBPSw5.p <- get_covariate_balance(pm.CBPSw5.p$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

mtext("Before\nRefinement", line = -1.8, at = 0.09, outer = TRUE, cex = 1)
mtext("Mahalanobis", line = -0.2, at = 0.26, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.42, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.59, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.75, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.92, outer = TRUE, cex = 1)

mtext("Standarized Mean\nDifferences for Protest", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# VSP
################################################################

pdf(file = "Figure9.1-V.pdf",   
    width = 13, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 6))


coba_pm.none5.v <- get_covariate_balance(pm.maha5b.v$att,
                                         data = data_final,
                                         use.equal.weights = TRUE,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.maha5b.v <- get_covariate_balance(pm.maha5b.v$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

coba_pm.ps.m5.v <- get_covariate_balance(pm.ps.m5.v$att,
                                         data = data_final,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.ps.w5.v <- get_covariate_balance(pm.ps.w5.v$att,
                                         data = data_final,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.CBPSm5.v <- get_covariate_balance(pm.CBPSm5.v$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

coba_pm.CBPSw5.v <- get_covariate_balance(pm.CBPSw5.v$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

mtext("Before\nRefinement", line = -1.8, at = 0.09, outer = TRUE, cex = 1)
mtext("Mahalanobis", line = -0.2, at = 0.26, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.42, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.59, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.75, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.92, outer = TRUE, cex = 1)

mtext("Standarized Mean\nDifferences for VSP", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()


################################################################
# Sanctuary
################################################################

pdf(file = "Figure9.1-S.pdf",   
    width = 13, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 6))


coba_pm.none5.s <- get_covariate_balance(pm.maha5b.s$att,
                                         data = data_final,
                                         use.equal.weights = TRUE,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.maha5b.s <- get_covariate_balance(pm.maha5b.s$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

coba_pm.ps.m5.s <- get_covariate_balance(pm.ps.m5.s$att,
                                         data = data_final,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.ps.w5.s <- get_covariate_balance(pm.ps.w5.s$att,
                                         data = data_final,
                                         covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                        "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                        "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                         plot = TRUE,
                                         legend = F,
                                         ylim = c(-1, 1))

coba_pm.CBPSm5.s <- get_covariate_balance(pm.CBPSm5.s$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

coba_pm.CBPSw5.s <- get_covariate_balance(pm.CBPSw5.s$att,
                                          data = data_final,
                                          covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                         "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                         "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                          plot = TRUE,
                                          legend = F,
                                          ylim = c(-1, 1))

mtext("Before\nRefinement", line = -1.8, at = 0.09, outer = TRUE, cex = 1)
mtext("Mahalanobis", line = -0.2, at = 0.26, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.42, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.59, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.75, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.92, outer = TRUE, cex = 1)

mtext("Standarized Mean\nDifferences for Sanctuary", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

##############################################################################
### COBA 2
##############################################################################

################################################################
# Protest 2
################################################################

pdf(file = "Figure9.2-P.pdf",   
    width = 13, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))


bscat_pm.5.p_maha.b <- balance_scatter(non_refined_set = pm.none5.p$att,
                                       matched_set_list = list(pm.maha5b.p$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

bscat_pm.5.p_ps_m <- balance_scatter(non_refined_set = pm.none5.p$att,
                                     matched_set_list = list(pm.ps.m5.p$att),
                                     data = data_final,
                                     covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                    "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                    "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                     main = "")

bscat_pm.5.p_ps_w <- balance_scatter(non_refined_set = pm.none5.p$att,
                                     matched_set_list = list(pm.ps.w5.p$att),
                                     data = data_final,
                                     covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                    "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                    "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                     main = "")

bscat_pm.5.p_cbps_m <- balance_scatter(non_refined_set = pm.none5.p$att,
                                       matched_set_list = list(pm.CBPSm5.p$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

bscat_pm.5.p_cbps_w <- balance_scatter(non_refined_set = pm.none5.p$att,
                                       matched_set_list = list(pm.CBPSw5.p$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.2.no_p_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Standarized Mean\nDifferences After Refinement", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Standarized Mean Differences Before Refinement", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# VSP 2
################################################################

pdf(file = "Figure9.2-V.pdf",   
    width = 13, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))


bscat_pm.5.v_maha.b <- balance_scatter(non_refined_set = pm.none5.v$att,
                                       matched_set_list = list(pm.maha5b.v$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

bscat_pm.5.v_ps_m <- balance_scatter(non_refined_set = pm.none5.v$att,
                                     matched_set_list = list(pm.ps.m5.v$att),
                                     data = data_final,
                                     covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                    "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                    "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                     main = "")

bscat_pm.5.v_ps_w <- balance_scatter(non_refined_set = pm.none5.v$att,
                                     matched_set_list = list(pm.ps.w5.v$att),
                                     data = data_final,
                                     covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                    "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                    "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                     main = "")

bscat_pm.5.v_cbps_m <- balance_scatter(non_refined_set = pm.none5.v$att,
                                       matched_set_list = list(pm.CBPSm5.v$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

bscat_pm.5.v_cbps_w <- balance_scatter(non_refined_set = pm.none5.v$att,
                                       matched_set_list = list(pm.CBPSw5.v$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.4b.no_v_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Standarized Mean\nDifferences After Refinement", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Standarized Mean Differences Before Refinement", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Sanctuary 2
################################################################

pdf(file = "Figure9.2-S.pdf", 
    width = 13, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))


bscat_pm.5.s_maha.b <- balance_scatter(non_refined_set = pm.none5.s$att,
                                       matched_set_list = list(pm.maha5b.s$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

bscat_pm.5.s_ps_m <- balance_scatter(non_refined_set = pm.none5.s$att,
                                     matched_set_list = list(pm.ps.m5.s$att),
                                     data = data_final,
                                     covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                    "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                    "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                     main = "")

bscat_pm.5.s_ps_w <- balance_scatter(non_refined_set = pm.none5.s$att,
                                     matched_set_list = list(pm.ps.w5.s$att),
                                     data = data_final,
                                     covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                    "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                    "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                     main = "")

bscat_pm.5.s_cbps_m <- balance_scatter(non_refined_set = pm.none5.s$att,
                                       matched_set_list = list(pm.CBPSm5.s$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

bscat_pm.5.s_cbps_w <- balance_scatter(non_refined_set = pm.none5.s$att,
                                       matched_set_list = list(pm.CBPSw5.s$att),
                                       data = data_final,
                                       covariates = c("cv1.lpop", "cv2.rur_p", "cv3c.c_lsh1", "cv4.nbi", 
                                                      "cv5c.1.bd_o", "cv6b.1.ce_o_r", "cv7b.3.pgf_cd_o_r", 
                                                      "cv8.3.no_s_str", "dv4b.2.1.farc_ckill_up1_r"),
                                       main = "")

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Standarized Mean\nDifferences After Refinement", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Standarized Mean Differences Before Refinement", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# End of this part
################################################################
